gr%% AudioPlayground
%  Script by adqeor@XJTU

%% 强制限幅引发失真

close all;

tSpan = 0.5;
Fs = 48e3;
t = 0:1/Fs:tSpan;
f = 440;

SC = 1;
clipRatio = .19;

y = SC*sin(2*pi*f*t);

yClipped = y;
yClipped( yClipped > SC*clipRatio) = SC*clipRatio;
yClipped( yClipped < -SC*clipRatio) = -SC*clipRatio;

figure; ax1 = subplot(211); ax2 = subplot(212);
plot(ax1, t(1:1000), y(1:1000)); ylim([-SC SC]);
plot(ax2, t(1:1000), yClipped(1:1000)); ylim([-SC SC]);

% soundsc(y, Fs); % soundsc 将音频信号摆幅缩放到正负 1 后播放
% pause(t_span); % sound 和 soundsc 是非阻塞播放, 接连两个 sound 会让两段声音叠起来
% pause(t_span);
% soundsc(yClipped, Fs);


figure; pspectrum([y; yClipped]', Fs, 'power', 'FrequencyLimits', [20, 8e3] ); % MATLAB 以列向量为信号通道划分 
legend({'y','yClipped'});

% figure('Name','y'); thd(yClipped, Fs); figure('Name','yClipped'); thd(yClipped, Fs);


%% 扫频测试
% 可用于测试室内混响
Fs = 40e3;
tSpan = 4;
t = 0:1/Fs:tSpan;

freq1 = 300;
freq2 = 400;

% 线性扫频, f1 扫到 f2
y = chirp(t, freq1, tSpan, freq2, 'logarithmic');

soundsc(y, Fs);

%{
% 这个也可
Fs = 40e3;
tSpan = 0.5;
t = 0:1/Fs:tSpan;

f0 = 200;
f1 = 400;
n = 20;

freqs = logspace(log10(f0), log10(f1), n);
% freqs = 12e3:.5e3:20e3;

for freq = freqs
	disp(freq);
	y = sin(2*pi*freq .* t);
	audioplayer(y, Fs).playblocking();
end
%}

%% 欧拉盘
% 可见 AV26257372

rawData = importdata('D:/AV26257372.mp3');
[p, f, t] = pspectrum(rawData.data(:,1), rawData.fs);
%
plot(p(:,1));
yticklabels(f);

%%

Fs = 40e3;

f_carrier = 5e3;
t_carrier = 1/f_carrier; % = doubled pulse positive span



%% Musical_Instrument.m
% 16 Mar, 2021 合并自 AudioPlayground

pianoTone =[	1		2		3		4		5;
						1		0.1		0.3		0.1		0.2
						0		0		0		0		0];

Fs = 40e3;

f0 = 440;
scale = [0 1]; % 向上下的八度

f = [];
for idx = 1:length(scale)-1
	f = [f, logspace(log10(f0*power(2, scale(idx))), log10(f0*power(2, scale(idx+1))), 8)];
end

tSpan = 0.2;
t = 0:1/Fs:tSpan;

y = [];
for idx = 1:length(f)
	y = [y, harmonics(f(idx).*pianoTone(1,:), pianoTone(2,:), pianoTone(2,:), t)]; % flip( gredge(length(t), 'ecl_4', length(t)*.4) )
end

sound(y, Fs);

%% Shepard Tones
Fs = 40e3;

f0 = 440;
scale = [-1 0]; % 向上下的八度

f = [];
for idx = 1:length(scale)-1
	f = [f, logspace(log10(f0*power(2, scale(idx))), log10(f0*power(2, scale(idx+1))), 8)];
end

tSpan = 0.2;
t = 0:1/Fs:tSpan;

y = [];
pianoTone =[	1		2		3		4		5;
						1		0.5		0.1		0.1		0.2
						0		0		0		0		0];
toneIdx = 1;
for idx = 1:12
	
	octaveIdx = 0;
	y_temp = zeros(size(t));
	while octaveIdx * 7 + toneIdx <= length(f)
		y_temp = y_temp + harmonics(f(octaveIdx * 7 + toneIdx).*pianoTone(1,:), pianoTone(2,:), pianoTone(2,:), t);
		octaveIdx = octaveIdx + 1;
	end
	
	y = [y, y_temp]; % flip( gredge(length(t), 'ecl_4', length(t)*.4) )
	
	if mod(idx, 2) == 0
		y = [y, zeros(size(t)), zeros(size(t))];
	end
	
	toneIdx = toneIdx + 1;
	if toneIdx == 8, toneIdx = 1; end
end

sound(y, Fs);